Microcanonical Approach to the Simulation of First-Order Phase Transitions. 
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A generalization of the microcanonical ensemble suggests a simple strategy for the simulation 
of first order phase transitions. At variance with flat-histogram methods, there is no iterative 
parameters optimization, nor long waits for tunneling between the ordered and the disordered phases. 
We test the method in the standard benchmark: the Q-states Potts model ((5 = 10 in 2 dimensions 
and Q = 4 in 3 dimensions) , where we develop a cluster algorithm. We obtain accurate results for 
systems with 10^ spins, outperforming flat-histogram methods that handle up to lO** spins. 
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Phase transitions are ubiquitous (formation of quark- 
gluon plasmas, evaporation/crystallization of ordinary 
liquids, Cosmic Inflation, etc.). Most of them are of 
(Ehrenfest) first order Monte Carlo simulations Q 
are crucial for their investigation, but difficulties arise for 
large system linear size, L (or space dimension, D). The 
intrinsic problem is that, at a first order phase transition, 
two (or more) phase coexist. The simulated system tun- 
nels between pure phases by building an interface of size 
L. The free-energy cost of such a mixed configuration is 
EL^~^ {E: surface tension), the interface is built with 
probability exp[— Z'L'^^^] and the natural time scale for 
the simulation grows with L as exp [rL^-i]. This disas- 
ter is called exponential critical slowing down (ECSD). 

No cure is known for ECSD in canonical simulations 
(cluster methods [1,|3| do not help), which motivated the 
invention of the multicanonical ensemble . The multi- 
canonical probability for the energy density is constant, 
at least in the energy gap e° < e < e'^ (e° and e'*: energy 
densities of the coexisting low-temperature ordered phase 
and high-temperature disorderedphase), hence the name 
fiat- histogram methods [^, 0, 0, S] ■ The canonical prob- 
ability minimum in the energy gap (cx exp[—SL^~^]) is 
filled by means of an iterative parameter optimization. 

In fiat-histogram methods the system performs an en- 
ergy random walk in the energy gap. The elementary 
step being of order L^^ (a single spin-flip), one naively 
expects a tunneling time from e° to e"^ of order L^^ spin- 
flips. But the (one-dimensional) energy random walk 
is not Markovian, and these methods suffer ECSD [lo| . 
In fact, for the standard benchmark (the Q = 10 Potts 
model Min D = 2), the barrier of 10'' spins was reached 
in 1992 ^, while the largest simulated system (to our 
knowledge) had 4 x 10'* spins 0. 

ECSD in flat histogram simulations is probably under- 
stood [13]: on its way from e'^ to e°, the system under- 
goes several (four in D = 2) "transitions" . First comes 



the condensation transition [lO|, lll| , at a distance of or- 
der from e'^, where a macroscopic droplet of 
the ordered phase is nucleated. Decreasing e, the droplet 
grows to the point that, for periodic boundary conditions, 



it reduces its surface energy by becoming a strip [l3| , see 
Fig. P (in D = 3, the droplet becomes a cylinder, then a 
slab [13|). At lower e the strip becomes a droplet of disor- 
dered phase. Finally, at the condensation transition close 
to e° , we encounter the homogeneous ordered phase. 

Here we present a method to simulate first order tran- 
sitions without iterative parameter optimization nor en- 
ergy random walk. We extend the configuration space 
as in Hybrid Monte Carlo [l3|: to our N variables, at 
(named spins here, but they could be atomic positions) 
we add N real momenta, Pi. The microcanonical en- 
semble for the {cTijPi} offers two advantages. First, mi- 
crocanonical simulations [isj are feasible at any value 
of e within the gap. Second, we obtain Fluctuation- 
Dissipation Eqs. dlHl]) where the (inverse) temperature 
f3, a function of e and the spins, plays a role dual to that 
of e in the canonical ensemble. The e dependence of the 
mean value (/3)e, interpolated from a grid as it is almost 
constant over the gap, characterizes the transition. We 
test the method in the Q-states Potts model, for which 
we develop a cluster algorithm. We handle systems with 
10^ spins for Q = 10 in D = 2 and for Q = 4 in I? = 3 
(where multibondic simulations handle N = 10"^ [l3|)- 

Let U be the spin Hamiltonian. Our total energy is 



N 



U. 



{e = £/N, u = U/N). (1) 



In the canonical ensemble , the {pi \ are a trivial gaussian 
bath decoupled from the spins. Note that, at inverse 
temperature (3, one has (e)^ = {u)^ + 1/(2/3) . 

Microcanonically, the entropy density, s(e, N), is given 
by (X]{(T }• summation over spin configurations) 

exp[Ns{e, N)] = J] E '^(^^ " ^) ' (2) 
or, integrating out the {pi} using Dirac's delta function, 
exp[iVs(e,7V)] = -§^5^^ J2 {e-u)^'^0{e-u) . (3) 
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The Heaviside step function, 9{e—u), enforces e > u. The 
microcanonical average at fixed e of a generic function of 
e and the spins, 0(e, {ci}), is (see Eq. ^ and 



(4) 



The Metropolis simulation of Eq. is straightforward. 
Calculating ds /de from Eq. ^ we learn that 31 1 



ds(e,iV) 
d^ 

/3(e;{aJ) = 



(/3(e;{aJ)), 

iV-2 
2iV(e - u) ' 



Fluctuation-Dissipation follows by derivating Eq. ([4]): 
d(0). 



de 



/do 
\ 9^ 



(0/3)e-(0)e(/5>. 



(5) 
(6) 

(7) 



As in the canonical case [l^, an integral version of ([7]) 
allows to extrapolate (0}e' from simulations at e > e': 



0{e';{a,})9{e' -u) 



e{e' - u) 



(8) 



For e < e', configurations with e < u < e' , suppressed by 
a factor (e'— are ignored in ([5]). Since we are lim- 
ited in practice to |e — e'| < yj{v?)e — {u)1/\d{u)elde\ ~ 
N~^/'^, the restriction e > e' can be dropped, as it is 
numerically negligible. 

The canonical probability density for e, Pp^\e) oc 
exp[iV(s(e,7V) - /3e)] follows from 0)^. 



logP^^'(e2)-logPr(ei)=^ 



de 



(9) 

In the thermodynamically stable region (i.e. d(/3)e/de < 
0), there is a single root of (/3)e — (3, at the maximum 

of -Pg^^ But, see Fig. [U in the energy gap {j3)e has a 
maximum and a minimum (L-dependent spinodals [ij), 
and there are several roots of {(3)e = (3. The rightmost 
(leftmost) root is e^(/3) (e5,(/3)), a local maximum of P^^^ 
corresponding to the disordered (ordered) phase. We de- 
fine e*i^{P) as the second rightmost root of (/3)e = (3- 
At the finite-system (inverse) critical temperature, (3^ , 



one has lli 



(L) 



which is 



equivalent, Eq. ^ and 20], to Maxwell's construction: 



= 



de (/3)e-/3e' 



(10) 



(for large N, (3^ - [3^ (xl/N 21]). Actually, at fixed 
e in the gap, also {i3)e tends to (3"^ for large N. In the 
strip phase it converges faster than see Table HI 
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FIG. 1: (Color online) Excess of (/3)e over /J^-^"" vs. e, for the 
Q = 10, D — 2 Potts model and several system sizes. Bottom: 
magnification for L > 512. The flat central region is the strip 
phase (the strip width varies at fixed surface free-energy). 
Lines (shown for L — 1024) are the two interpolations used 
for L > 512. We connect 3 independent cubic splines, in 
the strip pheise and in its sides, either by a linear function 
or by a step-like 1/100 power. Differences among the two 
interpolations are used to estimate the error induced by the 
uncertainty in the location of the strip-droplet transitions. 



In a cubic box the surface tension is estimated as |3S 
N 



de 



{0)e - Pi 



(11) 



L — > oo extrapolations 17°° — oc 1/L [23] are popular. 

As for the specific heat, for ^ oo the inverse func- 
tion of the canonical {e)p is the microcanonical {I3)e'- 



dp 



1 



1 



2(/3)2 d(/?)e/de 



CL{e). (12) 



For large N, ei(/3<P) , el{(3i^) , Ci(ei(/3<P)) , CL(ei(/3,^)) 
tend to e'^, e°, or the specific heat of the coexisting phases 
(we lack analytical hints about convergence rates). 

We now specialize to the Potts model [9|. The spins 
(Ti = 0; 1, . . . , Q — 1 , live in the N ~ nodes of a (hy- 
per)cubic lattice of side L with periodic boundary condi- 
tions, and interaction (< ij >: lattice nearest-neighbors) 



U 



<ij> 



(13) 



A cluster method is feasible. Let k be a tunable pa- 
rameter and w{e, u, k) = (e — u)^/'^~^ exp[KNu\9{e — u). 
Our weight is w{e, u, K)exp[— k[/] , see (|4]), or, introducing 
bond occupation variables, = 0, 1, and p = 1 — exp[K], 



w(e,u,K) Y\ [(1 -V)5n 

<i,J> 



(14) 
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FIG. 2: (Color online) L — 1024 equilibrium configurations 
for the ferromagnetic Q = 10 ,D = 2 Potts model with peri- 
odic boundary conditions, at the 2 sides of the droplet-strip 
transition, namely e = —0.809 (left) and e = —0.8 (right). 



which is the canonical statistical weight at = k [2J], 



but for the {riij} independent factor w{e,u,K). Hence, 
clusters are traced in the standard way, but we ac- 
cept a single-cluster flip Q with Metropolis probability 



p{e,K) = imii{l, w{e,u ,K)/w{e,u 



.initial 



t)}. Eqs. du- 



ll]) suggest that K = {(3)e maximizes p(e, k) (a short 
Metropolis run provides a first k estimate). We obtain 
(p(e, K))e > 0.99 for e < e'^, and still (p(e, K))e=e° > 0.78. 

We simulated the (Q = 10,1? = 2) Potts model [13], 
for L = 32,64,128,256,512 and 1024, sampling 0)e at 
30 points evenly distributed in —1.41666 < e < —0.45. 
For L — 512, we made 15 extra simulations to re- 
solve the narrow spinodal peaks (26 extra points for 
L = 1024). Our Elementary Monte Carlo Step (EMCS) 
was: max{10, Af/((A/')e(p(e, K))e)} cluster-flip attempts 
{J\f: number of spins in the traced cluster; it is of order 
one at e'^ and of order N at e°). So, every EMCS we 
flip at least N spins. For each e, we performed 2 x 10® 
EMCS, dropping the first 10% for thermalization. A sim- 
ilar computation was carried out for the {Q = A,D = 3) 
Potts model [11] (for details see Table |T] and 

Our (/3)e in Z? = 2 is shown in Fig. [TJ Data reweigth- 
ing ([8|) was used only to reconstruct the narrow spinodal 
peaks. To find the roots of {$)e = P, or to calculate 
the integrals in Eqs. (|10llip . we interpolated (/3)e using 
a cubic spline [33]. For L > 512 the strip-droplet tran- 
sitions produce two "jumps" in (/9)e, causing oscillations 
in the interpolation (Gibbs phenomenon), cured by either 
of two interpolation schemes, see Fig. [TJ 

We obtain I7^ eUf}^) , ei{(3^) , CL{e°M)) 

and CLie'liPc)) from the interpolation of (/3)e, and of 
d(/3)p/de, see ([7]). Statistical errors are Jack-Knife's 
[26j (the z-th block is obtained interpolating the i-th 
Jack-Knife blocks for {(3)e)- There are also interpola- 
tion and integration errors. Fortunately, errors of order 
e in e1{(3^) or e^(/3f') yield errors of order in (3^: the 
main error in is the quadrature error for (/3)e divided 



by the latent heat. On the other hand, e^(/3^) is near to 
the droplet-strip transition, and an error on it does have 
an impact on S^. 

In Table U are our results for (D — 2,Q — 10) and the 
known large L limits. A fit for c in P^-P^ = c/L^ ^ is 
unacceptable for L > 32 (x^/d.o.f. = 14.32/4), but good 
for L > 64 (x^/d.o.f. — 1.77/3): our accuracy allows to 
detect subleading corrections. A fit e'^(/3^) — e° — bi/L^ 
works only for L > 256 (xVd.o.f. = 1.90/2; for efiP^) 
we get x^/d.o.f. = 1.41/2). Rowevev, P^^"^'^ (see caption 
to Table m is compatible with P^ for L > 256. Then, 
the simplest strategy to get P^ and the latent heat is: 
(1) for L large enough to display a strip phase, locate it 
with short runs, (2) get /J*'''''?'^ accurately, and (3) find 
the leftmost (rightmost) root for (/3)e = fj'^^^P'^ . 

As for S^, the inequahty < 0.0473505 [13] (equal- 
ity under the hypothesis of complete wetting) was vio- 
lated by 1/L extrapolations performed with L < 100 
The reader may check (Table|I| that our data for L < 256 
extrapolate above 0.0473505, but drop below for L > 512. 
Indeed, the consistency of our results for P^ imply that 
the integration error for (/3)e is (at most) 2 x 10^® for 
L = 1024. Hence, the integration error for Sl is at 
most 10^'^. Adding it to the difference between the lin- 
ear and the step-like interpolation. Fig. [1] we obtain 
^L=i024 ^ 0.043(2), which is sHghtly below 0.0473505. 

As for (Q = 4, = 3), see Table U (3^^"p-L has con- 
verged (within accuracy) for L > 64. Hence, our pre- 
ferred estimate is P^ =0.6286206(10), that may be com- 
pared with Janke and Kapler's P^ = 0.62863(2) [l^. 
Accordingly, we find e°(/5'''"P"^) = -1.10537(4), 
gd(-^strip,L) ^ -0.52291(2), CL(e°(^"*"P'-^)) = 35.4(9), 
and CL(e'^(/3"*"P'-^)) = 4.24(18). The reader wiU note that 
P^^'^^ is far too high (for instance, from the x^/d.o.f . of 
the extrapolation P^ = P^ + cL^^). Therefore, the in- 
tegration error is ~ 4 x 10~® (larger than the statistical 
one) , which provides a bound for the error in the surface 
tension: Z'^=i28 = 0.0118(4). This is compatible with 
S^^^^, and provides a reasonable . 

We propose a microcanonical strategy for the Monte 
Carlo simulation of first-order phase transitions. The 
method is demonstrated in the standard benchmarks: the 
Q = 10, D = 2 Potts model (where we compare with exact 
results), and the (3 = 4, D = i Potts model. For both, 
we obtain accurate results in systems with more than 10® 
spins (preexisting methods handle 10^ spins). Envisaged 
applications include first-order transitions with quenched 
disorder [li', , colloid crystallization ^2^ , peptide ag- 
gregation ^] and the condensation transition [llj . 

We thank for discussions L. A. Fernandez (who also 
helped with figures and C code), L. G. Macdowell, W. 
Janke, G. Parisi and P. Verrocchio, as well as BIFI and 
the RTN3 collaboration for computer time. We were 
partly supported by BSCH— UCM and by MEC (Spain) 
through contracts BFM2003-08532, FIS2004-05073. 
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TABLE I: System size dependent estimates of the quantities characterizing the first order transition, as obtained for the 
Q = 10,-D = 2 Potts model (top) and (5 = 4,1) = 3 (bottom). Errors are Jack-Knife's. Also shown is = (/3>e=-o.95 (for 

D — 2) or ^sti'iP'L — (/3)^^_Q 7g4443 (for Z) = 3), in the strip phase. The oo^ row contains exact results [S^l and an inequality [27| ]. 
for D = 2, Q — 10. The results with superscript A{B) were obtained with the linear(step-like) interpolation scheme, see Fig. [l] 
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